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Abstract 

Multilevel Monte Carlo simulations of a BSCCO system are carried out including both Josephson 
as well as electromagnetic couplings for a range of anisotropics. A first order melting transition of 
the flux lattice is seen on increasing the temperature and/or the magnetic field. The phase diagram 
for BSCCO is obtained for different values of the anisotropy parameter 7. The best fit to the 
experimental results of D. Majer et al. [Phys. Rev. Lett. 75, 1166 (1995)] is obtained for 7 ~ 250 
provided one assumes a temperature dependence A^(0)/A^(r) = 1 — t of the penetration depth with 
t = T/Tc- Assuming a dependence A^(0)/A^(T) = 1 — t^ the best fit is obtained for 7 450. For 
finite anisotropy the data is shown to collapse on a straight line when plotted in dimensionless units 
which shows that the melting transition can be satisfied with a single Lindemann parameter whose 
value is about 0.3. A different scaling applies to the 7 = 00 case. The energy jump is measured 
across the transition and for large values of 7 it is found to increase with increasing anisotropy 
and to decrease with increasing magnetic field. For infinite anisotropy we see a 2D behavior of 
flux droplets with a transition taking place at a temperature independent of the magnetic field. 
We also show that for smaller values of anisotropy it is reasonable to replace the electromagnetic 
coupling with an in-plane interaction represented by a Bessel function of the second kind {Kq), 
thus justifying our claim in a previous paper. 
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I. INTRODUCTION 



High-temperature superconductors are materials of type II, that allow for partial magnetic 
flux penetration if the external field satisfies H^i < H < iJc2^'^'^- The flux penetrates the 
sample in the form of flux-lines (FL's), each containing a quantum unit 0o = hc/2e of 
flux. At low temperature the FL's form an ordered hexagonal lattice (Abrikosov lattice) 



due to their their mutual repulsion. The lattice constant is given by = y2(j)o/\/3B, 
where B is the magnetic field flux density. At high temperature and/or magnetic field this 
lattice melts due to thermal fluctuations^'^'^'^. Two of the most commonly used high- 
temperature superconductors are YBa2Cu307_5 (YBCO) and Bi2Sr2CaCu208+5 (BSCCO). 
They are very different in one very important respect - their anisotropy parameter 7, defined 
as 7^ = mz/m±, where rriz and m± denote the effective masses of electrons moving along the 
c axis and the ab plane, respectively. While for YBCO the anisotropy is somewhere between 
5-7, for BSCCO it is estimated to be between 10 to a 100 times larger, and estimates vary 
across the literature. What contributes to the discrepancy of different estimates is the fact 
that the anisotropy varies depending on the degree of doping. There is also the question of 
how to extract the value of anisotropy correctly from the exerimental results. Blatter et 
cite a range of anisotropics of 50-200 for BSCCO. References j^ls 10 1 cite values in the range 



of 140-160. However more recent acurate measurement the components of the London 
penetration depths Ac and Xab, the ratio of which is 7, are consistent with anisotropy in the 
range of 300-500 for optimally doped samples. The previously reported lower values may 
belong to overdoped samples or constituted only a lower bound. Anisotropy controls the 
amount of "wiggling" of a flux- line from plane to plane. In YBCO the FL's are more rigid 
while in BSCCO they are so loose that they are customarily referred to as a stack of two 
dimensional pancakes^^ (or droplet vortices) rather than FL's. 

Interaction between two FL's in YBCO is non-local. It is a screened Bio-Savart type 
of interaction where each segment of a FL interacts with every other segment of the same 
FL and all the other FL's. For segments oriented in the same direction the interaction is 
repulsive^: 

Here Sj(z) denotes the position of the i'th FL at elevation z along the 2;-axis, Eq = 0o/ (47rA)^ 
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is the line energy and A is the screening length (penetration depth). Considering two given 
FL's, it turns out though that it is a good approximation to replace the non-local interaction 
of a given line segment of one FL, with all segments of the other PL with a single interaction 
among segments belonging to the same plane. This interaction is given approximately hy^ 
2 Eq d KQ{Rij / X) where Rij is the distance between the two segments in the same plane and 
d is the thickness of each layer. The approximation is valid when the FL's do not deviate 
too much from straight lines which is a good approximation for YBCO in the "solid" vortex 
lattice phase, because FL's are stiff and do not wiggle too much. For each FL, there is also 
an elastic energy associated with its deviation from a straight line along the z-direction. 
The elastic energy of a flux-line in YBCO is approximately given by 



assuming the external magnetic field is aligned along the z-direction. The elastic coefficient 
(line tension) ei is equal to Eq ln(A/^)/7^ where ^ is the coherence length, and 7 is the 
anisotropy. In the discrete case this self-energy transforms into an attractive quadratic 
interaction between segments in adjacent planes. In this form the problem is equivalent to 
a system of bosons with repulsive interaction o^^i^^ . The term described in the last equation 
corresponds to the kinetic energy of the bosons which repel each other with a screened 
Coulomb interaction. 

For BSCCO the situation is different because each FL is represented more faithfully by a 
collection of pancakes. Each pancake interacts with every other pancake, but the interaction 
is different from the interaction among FL segments discussed above. The interaction can 
be shown to consist of two parts. The first part is called the electromagnetic interaction (or 
simply magnetic) and it exists even in the case that the layers of the materials are completely 
decoupled, so no current can fiow along the c-axis of the sample. A pancake vortex located 
in one plane gives rise to screening currents in the same plane as well as in all other planes. 
A second pancake vortex, located elsewhere, interacts with the screening currents induced 
by the first pancake^^. This interaction has been calculated by Clem and othersi^. Two 
pancakes in the same plane interact with a repulsive interaction while pancakes in different 
planes attract one another. If one considers a single pancake vortex and an infinite set 
of pancakes a distance R away stacked along the z-axis, then the interaction still sums to 
2 Eod Kq{R/ X) in the limit when d/X goes to zero (see Appendix A). Clen>i^ proceeds to 




(1.2) 
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show that if one has a straight array of pancake vortices along the z-axis, and one pancake 
of the stack is displaced a distance R in the lateral direction then the magnetic energy of 
the configuration increases by an amount 



where C is Euler's constant (=0.5772...). For large i? (i? ^ A), the modified Bessel function 
Ko decays exponentially and thus the energy increases like ln(i?/A). For small R the Bessel 
function can be expanded in a power series in R/ X 



and thus the electromagnetic energy behaves like R^ to leading order in R. 

The second part of the interaction among pancake vortices is the so-called Josephson 
interactio n^'^^1^^ . It results from the fact that there is a Josephson current flowing between 
two superconductors separated by an insulator and this current is proportional to the sine 
of the phase difference of the superconducting wave functions. The superconductors in the 
present case are the different Cu02 planes. When two pancakes belonging to the same 
stack and residing in adjacent planes move away from each other, the phase difference that 
originates causes a Josephson current to begin flowing between the planes. This results in 
an attractive interaction between pancakes that for distances small compared to Vg = 'yd is 
approximately quadratio^*i^ in the distance. When the two adjacent pancakes are separated 
by a distance larger than Vg, a "Josephson string" is formed, whose energy is proportional 
to its length^. 

When the anisotropy is not too large, the Josephson coupling among adjacent pancakes, 
which are loosely belonging to the same "flux-line", dominates over the electromagnetic 
interaction, and the later can be neglected. The ratio of the coefficients of the quadratic 
terms in the effective electromagnetic interaction (as mentioned above) and the Josephson 
interaction goes roughly like 7^((i/A)^ (where d/A ~ 1/120 for BSCCO at T = and even 
smaller at higher temperatures). Thus for anisotropy 7 = 50 we get a factor of 0.25 or less (a 
somewhat more precise estimatei^ gives a ratio of about 0.1). Thus the magnetic interaction 
is small compared to the Josephson interaction for anisotropies in the range of 7 = 50 — 100. 
For samples with 7 = 200 these interactions are already comparable. For large values of R 
the magnetic interaction increases logarithmically and the Josephson interaction increases 





MR/X) 



ln(i?/2A)(l + i?V4A2 + ■■•)- C + R^il - C)/4:\^ + ■ ■ ■ , 



(1.4) 
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linearly so the electromagnetic interaction is always negligible. The key to the estimate 
given above is to consider not just two pancake vortices but a whole line with one displaced 
pancake. This argument is valid if the deviations of the vortices from straight lines are not 
too large. On the other hand, the electromagnetic interaction starts to be important for 
anisotropics which are significantly larger than \/d which for BSCCO is about 120. 

It is the aim of this paper to include both the electromagnetic interaction and the Joseph- 
son interaction among pancakes and to see what is the combined effect on the phase diagram 
of the melting transition and the energy jump across the transition. The electromagnetic 
interaction will be included fully in the sense that we will not make the approximation that 
the pancake stacks are nearly straight and hence the electromagnetic coupling will not be 
replaced by an in-plane effective coupling. Numerically, much of the past work on BSCCO 
has been confined to the X-Y mode l -^^i^^i^^i^^i^^i^^ and Bose modeU^^, both of them treat 
the electromagnetic coupling imprecisely by including it as an effective in-plane interaction. 
Recently in several papers using the Langevin simulation method^SiSLS^^ the electromagnetic 
coupling has been fully taken into account. Unfortunately, these papers completely neglect 
the Josephson coupling which can hardly be justified. Also some of these papers work with a 
small system size like 5 to 10 planes along the z-direction, and some do not even use periodic 
boundary conditions in the z-direction. In this paper we carry out Monte Carlo simulations 
of a BSCCO system consisting of 20 — 36 planes and periodic boundary conditions are used 
in all directions, including the z direction. Thus a pancake would interact with an infinite 
number of pancakes through the images under the periodic boundary conditions. For this 
an efficient way to sum over the interaction is required. We derive a formula for summing 
over the logarithmic interaction in Appendix B. 
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II. THE MODEL 



The starting point is the Lawrence-Doniacb^^ Gibbs free-energy functional, 
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(2.1) 



where ipn represents the superconducting order parameter in the n*^ Cu02 layer, a^^^ is 
the vector potential in the plane, and d is the thickness of the insulating layers, b is the 
local magnetic field and H the externally applied field. The usual 3D integration of the GL 
theory has been replaced by a summation over all the superconducting layers along with a 
2D integration over the superconducting planes. We set ipn = |^n|exp(z0„) and, working 
in the London approximation, we drop the term alV'nP + /3|'?/'n|^/2 because it gives only a 
constant contribution. Then we get 

g= I d'R- 
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Minimization with respect to (i.e. ^ = 0) and ay (i.e. |^ = 0) gives 



(2.2) 
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X'Aa^'^ =dT6,{z-nd) a(2) + |iv(2V„ , (2.4) 
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where A stands for the 3-dimensional Laplacian. Minimization with respect to gives 
End 2m 27r . ,^ 

where 
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is the gauge invariant phase difference between the layers n and n + 1. Eq. ()2.5|) imphes 

An 

Aa^ = — jjsin($„,„+i). (2.7) 
c 

where 

is the Josephson-couphng current density between layers. Minimization with respect to 
gives 

A(2)0„ + |^V(2) ■ a'^' = 4h [sm($„,„-i) - sin($„+i,„)] • (2-9) 

Eqs. ()2.4|) . ()2.7p and ()2.9|) are to be solved with the appropriate boundary conditions and 
the solution must be substituted back into Eq. ()2.2j) to obtain the Gibbs free-energy, and 
thus the strength of the interaction among pancake solutions. We also see that in the limit of 
infinite anisotropy the right- hand- side of equations ()2.7j) and ()2.9j) tend to zero. An isolated 
pancake residing in plane n is a singular solution of the equation for the phase of the wave 
function which satisfies 

V(^)0„(R) = _ "><(R-y , (2.10) 

where R is a two dimensional vector in the plane and Hn denotes the center of the pancake. 
By n we denote a unit vector in the z-direction. Thus as one fully encircles the pancake the 
phase (pn changes by 27r, and is singular at the center of the pancake. In the case when the 
Josephson coupling is totally neglected, i.e. for 7 — oo, the full solution of Eqs. ()2.4|) . ()2.7|) 
and ()2.9p can be found and from it one can easily obtain the magnetic-field in real space^^. 
It is not a trivial matter to switch from the variables 0„ and a to pancake variables and 
express the free energy in terms of the latter. This transformation can only be implemented 
approximately. In the following we first summarize the known results for the Josephson 
interaction and the electromagnetic interaction and then proceed to combine them together 
into a single algorithm. Most papers consider only one type of interaction, Josephson or 
electromagnetic in the limit that one dominates over the other. 

A. Josephson Coupling 



We keep the same Josephson coupling as in Ref. |25l This coupling is strongly dependent 



on the anisotropy parameter 7. Consider two adjacent pancakes, belonging to the same PL, 
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residing in the m and m + 1 planes respectively, such that their centers are displaced and 
not located on top of each other. Assume that the pancake on the m + 2'th plane is located 
on top of the m + 1-plane pancake, and the pancake residing in the m — 1 plane is at the 
same position as the m-plane pancake. Thus (pm+i 7^ 0m but (f)rn+2 = 4>m+i and (pm-i = 4>m- 
This assumption is made in order to trancate the infinite set of couples equations^, and in 
real situations may constitue an approximation. Denoting $m+i,m = <Pm+i — <Pm simply by $ 
we see that writing down Eq. ()2.9|) for the m and m + 1 planes respectively and subtracting 
one equation from the other, one obtains 

A(2)<I> = ^sin($), (2.11) 

where Vg = 'yd is the relevant screening length of the problem. Note that the screening 
term due to the vector potential has been neglected since it is negligible on length scales 
-R ~ ^ 7A. Equation 1)2.111) is the famous Sine Gordon Equation. Once its solution is 
obtained, it needs to be substituted in the Lawrence- Doniach Gibbs free-energy. This results 
in a contribution of the fornt^ii 

G.j = ^ ! d2R'[l - cos($(R')). (2.12) 

If we denote the separation between the two pancakes by R one can show that for R <^ 
R' <^ Tg the solution of Eq. ()2.11|) is given simply by 

$(R') = i?sin(^')/^', (2-13) 

since the right-hand-side can be neglected in this region. Here 9' is the azimuthal angle 
in the plane. For R' ^ Tg the sin($) can be replaced by $ and we see that the solution 
decays exponentially with a screening length Tg. Thus substituting Eq. ()2.13|) into Eq. ()2.12p . 
expanding the cosine to quadratic order and cutting off the integration at a large distance 
R' = Tg and small distance R' = R the interaction energy becomes 



dsQ , (rg\ f R 



SAR)-"-fl^[:i)[yJ, (2.14) 

SO we see that it is approximately proportional to R^ . On the other hand when the separation 
between the centers of the two pancakes becomes larger than Tg then a Josephson string is 
formed between the m'th and the m + I'th planes in a direction parallel to the planes. The 
energy of the Josephson string is proportional to its length which is equal to R, the pancakes' 
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separation. The calculation of this energy is rather involved and discussed by Clem, Coffey 
and Haoi^. The final result is 

Gj{R) = dso (^1.12 + In Q^^ ^. (2.15) 

The question now arises how to match these two interaction potentials valid for R <^ rg and 
R ^ Vg, one behaving quadratically in R and one linearly in R. One such extrapolation 
was given by Ryu, Doniach, Deutscher and Kapitulnik (RDDK)^^, who achieved a matching 
by keeping the coefficient of the linear term in R as given in Eq. ()2.15p in both regions, 
choosing the matching point to be at i? = 2rg, and subtracting a constant so that the 
two expressions would vanish at the matching point. Keeping the same constant in both 
expressions assures that the first derivative is continuous at the matching point. RDDK 
also replaced the constant 1.12 by 1. There are ways to improve the extrapolation^- but 
they do not change significantly the results of the simulations performed using the RDDK 
formula. Unfortunately there was a mistake by a factor of 7i/2 in the RDDK formula 
that we corrected below that has the effect of renormalizing the anisotropy parameter by 



about \prij2 ^ 1.25 since the major contribution to the simulations come from the region 
R < 2rg. To summarize, the London free-energy for inter-layer (IL) Josephson coupling is 
given approximately by: 



^>.L(R.™,R..^i) = ^(l + lnQ 
for |Ri,m - Ri,m+i| < 2rg , and 

S>/L(R*,m, Ri,m+l) = ^Q^2X2 ^ \ d 



(|Ri,m Rj,m+l|) 
^'9 
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(2.16) 



A 



i,m+l I 
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(2.17) 



for |Rj,m — Rj,m+i| > 2rg. Here the position of a pancake is specified in terms of cylindrical 
coordinates. Thus the position of the z'th pancake in the m'th plane is given by (Ri^m; md), 
where Ri,m is a f"wo dimensional vector in the ab plane. The index i labels the FL that the 
pancake is a part of. We have considered only pancakes belonging to the same FL. This is 
because for large separations by definition the Josephson string is formed among pancakes 
belonging to the same FL. We do check for the energy of all nearest neighboring pairs when 
deciding how to connect pancakes and we allow the process of flux "cutting" and switching. 
In the case that pancakes belonging to different FL's approach each other to a distance 
much smaller than Vg, than there will be an interaction of the form of Eq. ()2.14|) . however 



this occurrence is rather rare for magnetic fields of the strength considered here and hence 
neglected by RDDK and also in this work. 

When using Eqs. (j2.16|) and (j2.17|) it is necessary to specify the temperature dependence 
of A. Different choices for this dependence are found in the literature. In this work we used 
the same choice as in Refs.i^^i^iSi^ motivated by Ginzburg-Landau (GL): 

^ = l-r/T.. (2.18) 

Some author o^^i?^ use a dependence of the form 

^ = 1 - {T/T^r. (2.19) 

Recent experiments^ show that the temperature dependence of the london penetration 
depth is not universal and depend on the amount of doping and the sample history. In Fig. 

we show the temperature dependence of the three samples investigated in Ref . • Also 
shown on the figure is the two fiuid dependence^ A^(0)/A^(T) = 1 — (T/Tc)^ which normaly 
applies in the opposite limit of very small GL ratio k = A(0)/^(0), unlike high-Tc materials. 
Displayed also is the weak coupling, clean limit BCS curve. Both curves are lying above the 
experimental data. We added to the original figure the two behaviors given in Eqs. ()2.18p 
and ()2.19|) to show that the experimental data actually falls between theses two curves, 
so they give reasonable lower and upper bounds to the experimental results. In the result 
section we discuss how this choice of temperature dependence affects the comparison of the 
simulation results with experiments. 

B. Electromagnetic Coupling 

This coupling can be obtained from the Lawrence- Doniach model discussed at the begin- 
ning of this section by putting 7 = oo, which eliminates the Josephson coupling altogether. 
Extensive calculations can be found in the literatur o^i^i^^ . For the in-plane interaction be- 
tween two pancakes one finds, 

U(i?,„0) = 2deo + ^E,) , (2.20) 



2A / Rij 2A 



where Rij = |Ri,m — ^j,m\ is the radial distance in cylindrical coordinates and "m" denotes 
the index of the plane. 
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FIG. 1: Temperature dependence of the normalized penetration depth for 3 experimental samples 
as given in Ref. I35I . For comparison the two-fluid and BCS result for clean superconductors in the 
weak coupling limit are shown. We also added the linear and quadratic behavior as given in Eqs. 
1)2. 18(1 and ()2.19|) in the text, represented by dashed and dot-dashed lines respectively. 

The interaction between two pancakes situated at different planes (R-i^m, fnd) and 
{Rj,n, nd) is given by, 

U(i?,„z) = -^(^exp(-|^|/A)ln^-E2) , (2.21) 

where Rij = |Ri,m — Rj,n\, and z = {m — n)d. 
In the above equations we defined 

E^= dpexp{-p/\)/p, (2.22) 

J Rij 

POO 

E2= cipexp(-v/?Tp2/A)/p, (2.23) 

C is a constant of the order of the system's size that cancels out upon taking energy differ- 
ences. 

III. NOTES ON THE SIMULATIONS 



We work with M rhombically shaped cells stacked on top of each other in z direction. 
All of these cells have periodic boundary conditions in x and y directions. Each one of 
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these rhombic units cells also repeats itself every M plane due to the periodic boundary 
condition in z direction as well. Thus we have periodic boundary conditions in all directions. 
There are a total of pancakes in each of the M cells. We work with two system sizes 
to safeguard against any possibile finite size effects. While working with = 36 we chose 
M = 25 in most cases, except for 5 = 100 G and 7 < 150, where the number of planes was 
increased from 25 to 36. For low fields and anisotropics the entanglement length along the c 
axis becomes large, and hence in order to observe a sharp transition a larger system size in 
the z-direction is needed. Similarly, for 64 pancakes in each plane, we usually work with a 
total of 20, 25 or 36 planes depending on the values of parameter B and 7. In all the cases 
considered in this paper we always had a total of at least 900 beads (A^ = 36, M = 25) and 
a maximum of 2304 beads (A^ = 64, M = 36). A FL consists of one pancake from each and 
every plane. Pancakes belonging to a given FL were tracked with pointers and linked lists^. 

Pancakes were moved by either Metropolis algorithn>2i or its advance form, multilevel 
Monte Carlo (MMC)^. Which method to employ in a particular case usually depends upon 
the anisotropy and the magnetic field, as described below. 

For most anisotropics and magnetic fields, the MMC technique was used. In MMC we 
update several pancakes spread over many planes at once. Thus for the lower anisotropics 
7 = 125 and 150, and B = 100 G, a total of 15 to 20 beads spread over 5 planes were used 
to update the system. For other anisotropics and fields (7 > 250 and B < 900 G) it suffices 
to use just 3 planes in the MMC technique. 

For different parameter ranges, one needs to use different methods of updating the FL's. 
For example, one can not use MMC method for 7 = 00 since there is no natural choice to 
generate paths sampled with a free Gaussian distribution. Even at higher values of 7, such 
as 500, the MMC method using 5 planes would be slow and inefficient. Thus for 7 > 375 
and B > 300, we move only one pancake at a time using the Metropolis algorithm. Flux 
cutting was implemented by using a different kind of move to allow for the large wiggling 
of the FL's (and thus to avoid bias towards straighter FL configurations). In this move 
the two ends of neighboring lines were switched by considering only the relevant Josephson 
interactions without attempting to displace any of the pancakes involved. 

Just like the case with high anisotropics discussed above, in the MMC implementation 
also we allowed for large wiggling of FL's along the ^-direction by introducing the process of 
flux cutting. Starting from a few neighboring lines, we cut chunks of FL's spreading over a 
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number of planes. New paths between the starting and ending positions of the PL's thus cut 
were made using a random walk through the space of permutationsi^, and the subsequent 
use of the bisection method^. 

Of course there are ranges of parameters where more than one method of update can be 
employed (e.g. Metropolis method, MMC with 3 planes or MMC with 5 planes.). In these 
cases of overlap, relevant methods were found to lead to the same result, as they should. 
Further details on the simulation technique are supplied in Ref. |25. 

The logarithmic part of the in-plane and out of plane electromagnetic interaction (see 
Eqs. ()2.2()|1 and ()2.2Hl ) was handled by analytically summing over the interaction as shown 
in the Appendix B. The integrals Ei and E2 were evaluated numerically taking into account 
the periodic boundary conditions in all directions. This is accomplished by considering 
images of pancakes in all directions. 

In the simulations as T is raised it is B that is kept fixed, not H, as done in experiments. 
This causes the phase transition to be less sharp, especially at low fields (where the phase 
boundary is flatter), as is evident from Fig. [21 The figure shows both the H — T and B — T 
phase diagrams schematically. In the B — T phase diagram there is a region where both 
the vortex-solid and vortex-liquid phases coexist. The paths corresponding to increasing 
temperature at constant B (MC) and constant H (experiment) are shown. 

The broadening of the phase transition can be avoided by using 'isobaric" Monte Carlo 
simulations^^, but this was not done in the present work. 

We measured the following physical quantities. For details the reader is referred to our 
earlier work^. 



A. Energy 

An expression for energy can be obtained from 

E = kT'^Mm,(3,N)). (3.1) 

Due to the internal temperature dependence of A on T, one gets a very complicated expres- 
sion for energy (not written down here). However, a simplified expression for energy can 
be obtained under the assumption that < A, since in this case it is justified to ignore 
the internal temperature dependence of A on T while taking derivatives with respect to the 
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FIG. 2: Schematic phase diagram showing the vortex-soHd (VS) and votex-hquid (VL) phases in 
the H — T and B — T planes. The path the system traces as the temperature is raised at fixed B 
(MC) or fixed H (exp) are indicated 

temperature^^. The energy expression obtained for the case when ao < A is given in Ref. 
25l . This expression, however, was seen to work well even for the cases where Oq ~ A and 
the difference between the simplified and the exact energy calculations was found to be 
insignificant for all cases except for very low values of B such as 40 — 80 G (we used these 
values for only 7 = 00). We have used the simplified expression for the energy in all cases. It 
would not affect the melting transition in anyway. The only change will be that the energies 
obtained will be off by a few percents for the case when the magnetic field is very small. 



B. Translational structure factor 



The translational structure factor S'(Qi) is defined as. 



^(Q 



i) = ]i^(E^^^''"^'''''""''^""^7' (3.2) 
\ / 

where (...) stands for the MC average, and Qi stands for a reciprocal lattice vector corre- 
sponding to the first Bragg peak and is given by 

Qi = r^(ei - e2Cos6l), (3.3) 

ao sm t7 
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where 9 = 7r/3, is the nearest neighbor distance and ei^2 are the unit vectors along the 
hexagonal unit cell such that 

ei • 62 = cos 6'. (3.4) 



C. Line entanglement 

As we allow permutations of FL's, we can define a number N^/N as that fraction of 
the total number of FL's which belong to loops that are bigger than the size of a "simple" 
loop. A simple loop is defined as a set of M beads connected end to end, M being the total 
number of planes. Loops of size 2M, 3M... start proliferating at and above the transition 
temperature and in the corresponding 2D boson system this proliferation is related to the 
onset of the superfiuidity. 

Some of the other important parameters were taken as follows: Aq = 1700 A, = 15 A 
and Tc = 90 K. 

IV. RESULTS 

A. Josephson and electromagnetic coupling 

In this section we discuss the results when Josephson as well as electromagnetic couplings 
are included in the expression for the free-energy functional. Three different quantities, the 
translational structure factor at the first Bragg peak S'(Qi), the energy E and the line 
entanglement N^/N were monitored. Simulations were done for four different anisotropy 
parameters 7 = 125, 250, 375 and 500. In addition we carried out simulations for two different 
temperature dependence of the penetration depth. The results for the first temperature 
dependence of A, namely A^(0)/A^(T) = 1 — T/T^, are shown in Figs. Q-®- The results 
for the second dependence of A, i.e. A^(0)/A^(T) = 1 — (T/TcY will be discussed later, in 
the context of the phase diagram. 

To check against possible finite size effects we worked with two different system sizes, 
namely 36 and 64 FL's, as discussed in the previous section. The structure factor at the 
first Bragg peak, for the two different sizes is shown in panels (a) and (c) of Figs. OEl It is 
clear that the transition temperature is unaffected by the choice of the system size. Similar 
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FIG. 3: 7 = 125. For various fields S = 100 G (filled squares), 300 G (open circles), 500 G (filled 
circles), 700 G (open triangles), 900 G (filled triangles) and 5000 G (lower triangles), the following 
quantities are shown: (a) the translational structure factor at the first Bragg peak for = 64 
FL's (b) the translational structure factor at the first Bragg peak for = 36 FL's (c) Energy 
for = 64 FL's (energy is given up to an additive constant, which is not important.) (d) Line 
entanglement for = 64 FL's. 

agreement was seen in plots of energy vs. temperature and the graphs of entanglement vs. 
temperature. These later comparisions are not displayed. 

A first-order transition (POT) is seen for all anisotropies, which is in agreement with 
numerous experiment a l^ i ^ i ^i ^ as well as numerica l^^i^'^'i^^i^^i^^i^^i^^i^'^ studies on type-II super- 
conductors. The location of the POT is inferred from a sharp decay in S'(Qi) and a sharp 
rise in the line entanglement and a discontinuous jump in energy E. Except for the case 
with 7 = 125 and B = 100 G, the transition is sharp and can be easily located. In all the 
cases, we see a discontinuous jump in E. The size of the jump can be determined in the 
following way. Take two sets of points on ii^ vs. = T/(l — T/Tc) graph, one in the vortex 
solid phase and the other in the vortex liquid phase. This can be done by using 5'(Qi) or 
Nf,/N vs. Tr graphs. We fit a straight line to the first set of points and another straight line 
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FIG. 4: 7 = 250. The same quantities are shown as in Fig. |21but 7 = 250 here. 
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FIG. 5: 7 = 375. The same quantities are shown as in Fig. |31but 7 = 375 here. 
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FIG. 6: 7 = 500. The same quantities are shown as in Fig. |21but 7 = 500 here. 

to the other set of points. The two hnes are extended up to the transition temperature and 
the jump in E is read off. 

It is clear from the E vs. Tr graphs that for a given anisotropy one gets lower jumps in 
E at higher magnetic fields. This effect is more pronounced for higher anisotropics. 

B. Only electromagnetic coupling 

Many simulation studies completely neglect the Josephson coupling^^^ and work only 
with the electromagnetic coupling. In this section we discuss the results obtained with 
neglecting the Josephson interaction (7 = 00). The results are shown in Fig. [7| The 
results can be compared with a recent simulation study of the same system using substrate 
model^^^. The phase boundary obtained by Dodgson et al. falls almost on top of the phase 
boundary obtained in this paper. Another important feature of this study is the increasing 
energy jumps toward lower magnetic fields as well as almost diminishing jumps in E toward 
very high magnetic fields {B > 900 G). The entropy jump calculated using the formula 
As = AE/T shows an increasing trend towards higher temperature (not shown here) again 
in agreement with the references cited above. 
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FIG. 7: 7 = oo (no Josephson coupling), (a) the translational structure factor at the first Bragg 
peak for B = AO G (filled squares), 60 G (open circles), 70 G (filled circles), 80 G (open triangles), 
(b) Energy for the same fields as in (a) (c) the translational structure factor at the first Bragg peak 
for B = 100 G (filled squares), 300 G (open circles), 500 G (filled circles), 700 G (open triangles), 
900 G (filled triangles) and 5000 G (downward triangles), (d) Energy for the same fields as in (c). 

Also, for the very high magnetic fields, the transition takes place at T2D = 16 K, indepen- 
dent of the magnetic field. This result is in agreement with many theoretical studiea^ii^. In 
this regime of fields each set of pancakes in a plane melts irrespective of the other planes and 
shows a characteristic 2D melting transition at T2D ~ 16 — 19 K. This kind of 2D melting 
is expected at high fields as the in-plane pancakes get closer to each other and hence the 
in-plane interaction keeps getting stronger compared to out of plane interaction and the 
different layers start to melt independent of each othei"^. 

Comparing results with the previous section, it is clear that the Josephson interaction, 
present at finite anisotropy, has the effect of reducing the energy jumps compared to the 
case when only the electromagnetic coupling is included. 
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C. Phase diagram 



Based on the results from the previous two sections, we obtain a phase diagram for 
BSCCO which is shown in Fig. |H1 The diagram shows the melting curves for various 
anisotropics, including the case when there is no Josephson coupling. On the same dia- 
gram we also show two experimental lines, one for the melting transitior^*^ and the other 
for the irreversibility line^°. The experimental irreversibility line falls very close to the simu- 
lated 7 = 250 phase boundary. The actual experimental melting curve does not seem to fall 
on any of the melting lines obtained from the simulations. The experimental melting line 
bends downward toward lower temperatures and this feature could not be seen here. This 
may the result of point defects present even in pristine samples, an effect reported in recent 
simulationg|22i2i. 

In Fig. El we show the phase diagram obtained by assuming a different dependence of A 



on temperature, namely A(T) = A(0) /a/1 — T^/T^. For this choice the irreversibility line 
falls somewhat below the 7 = 500 simulated phase boundary, leading to an estimate of 
7 ^ 450. As discussed above the temperature dependence should be somewhere in between 
these dependences so an estimate of 7 = 300 — 400 for the experimental sample is probably 
reasonable. 

We have found that the data for the phase boundaries can be collapsed on a single curve 
for the case of finite anisotropy even when using different temperature dependences of A. 
When we plot the variable -87^ at the melting transition versus the variable kT / {eo{T)d) all 8 
curves fall on top of each other. Koshele\i2i argues that when Josephson coupling dominates 
over electromagnetic coupling the phase boudary is determined by a single dimensionless 
function of the dimensionless parameters kT/e^d where Sq{T) = 0Q/(47rA(T))^ and Tg/a^ oc 
572. In Fig.dTHlwe plot IniB'y^) versus ln(A;T/(eo'^)) and the data collapses to a straight 
line with slope (—2). This suggests that the transition is given by a single relation 



A 



(,V2eod) (ao) ^^'^^ 



as is known to be the case for YBCO^^. From a simple cage model (see e.g. Ref. 1361 ) 
it follows that the Ac = c|, where cl is the Lindemann coefficient. This means that the 
transition occurs when the mean square devaiation of the FL from a straight line exceeds 
clOo, where is the lattice constant. Since we find that Ac ~ 0.1, it appears that cl ~ 0.3 
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FIG. 8: Phase diagram obtained by using A(r) = A(0)/yi - T/Tc. Melting transitions for various 
anisotropies 7 = 125 (filled triangles), 7 = 250 (open circles), 375 (filled circles), 500 (open 
triangles), 00 (filled squares) and experimental(inverted triangles and diamonds) are shown. 

(It was estimated to be 0.25 for YBCO^^). Notice that actually the expression for A should 
be kT / {aQ^J2£leQ) and we have used Si = £0/7^- Had we included the factor \n{X/d) in 
El we would have obtained Ac = O.l/^/b ~ 0.045 and cl ~ 0.2. Thus cl turns out to be 
comparable to that found for YBCO. This phase diagram shows the importance of keeping 
the Josephson coupling. Even if the anisotropy parameter is as large as 7 = 500, it is still 
not appropriate to neglect the Josephson coupling entirely. 

For 7 = 00 there is another scaling which collapses the datai^^^ for the two kinds of 
temperature dependence we considered, namely plotting B/Bx at the melting transition vs. 
kT/{eod), where Bx = 0o/A^(T). This is an approximate scaling of the electromagnetic 
interaction valid when the dependence on the small parameter d/X can be neglected. In 
Fig. (fTTj) we show the collapsed data for two different temperature dependences and also 
compare with the result of Dodgson et a^^^. 
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FIG. 9: Phase diagram obtained by using A(r) = A(0)/a/1 - TV??: Melting transitions for 
various anisotropies 7 = 125 (filled triangles), 7 = 250 (hollow circles), 375 (filled circles), 500 
(hollow triangles), 00 (filled squares) and experimental melting transition (inverted triangles) and 
the experimental irreversiblity line (diamonds) are shown. 

D. Comparison: local and nonlocal interactions 

It is easy to show that the expression for the interaction of a single pancake with an infinite 
set of pancakes stacked on top of each other reduces to the modified Bessel function of second 
kind Kq{R/\) where R is the distance of the pancake from the stack of the pancakes (see 
Appendix A). At lower temperatures, magnetic field and anisotropy, the FL's are straight so 
it has been considered a good approximation to replace the nonlocal interaction between the 
pancakes with a local one of type mentioned above. For example we have used this model 
in a recent simulations^. In this section we compare the results obtained by replacing the 
nonlocal electromagnetic interaction with a local interaction of the Kq type. The results are 
shown in Fig. assuming the temperature dependence A^(0)/A^(T) = l — t. In Fig. ITW a) 
we show the phase lines for 7 = 125 and 250. The phase lines for 7 = 125 from the two 
models fall almost on top of each other. Even for 7 = 250 the deviation is still small. This is 
reasonable, as one expects that the FL deviates less from a straight line configuration when 
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FIG. 10: Data for finite anisotropies 7 = 125 ( triangles), 250 (circles), 375 (squares) and 500 
(inverted triangles) from Fig. |H1 and Fig. El collapse onto a single curve. Filled symbols are for the 
A2(T) = A2(0)/(l - T/Tc) dependence while hollow symbols are for X^{T) = A2(0)/(1 - T^/T^) 
dependence. Some spread in the data is attributed to the fact that this scaling is not valid at 
high values of anisotropies where Josephson coupling becomes comparable to the electromagnetic 
coupling. 

7 is small. In Fig. IT^T b-d) we compare the results obtained here with a previous paper— 
for 7 = 125 and B = 125 G. There is a small shift in transition temperature of about 3 K. 
It is clear that if one shifts the >S'(Qi) line obtained with Kq{R/X) interaction by around 3 
K to the right, then the two curves will fall almost on top of each other. 

Fig. IT^ b) compares the energy jumps for the same two case as in Fig. IT^ b). Except 
for slight shift in the transition temperature, the jumps in energy are identical. Fig. fT^ c) 
shows the line entanglement and as expected we again see a slight shift in the position where 
line entanglement takes place, consistent with Fig. ITW b) and lTW c). 

This justifies our claim that for the smaller anisotropies (7 ~ 125) it is a good approx- 
imation to replace the non local electromagnetic interaction with a local one of the type 
MR/X). 
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FIG. 11: Data for 7 = oo from Fig. |H1 (triangles) and Fig. IHl (squares) collapse onto a single curve. 
Also shown is the melting line obtained by Dodgson et al."^"^ (circles). 

V. CONCLUSIONS 

In this work we introduced a model which incorporates the electromagnetic interaction 
and Josephson interaction among pancake vortices in the highly anisotropic BSCCO in 
a more systematic way than previously done. Instead of approximating the effect of the 
electromagnetic interaction by an in-plane repulsive interaction between pancakes that is 
given by a Bessel function of the second kind we treat this interaction exactly and in addition 
we include the Josephson coupling. Treating the electromagnetic interaction as an effective 
Bessel function interaction is valid only if FL's do not deviate too much from straight lines. 
These distortions can be quite large close to the melting transition when the flux cutting 
mechanism starts to proliferate. Thus in the present model interaction among pancakes is 
taken in a more realistic manner. There are still some approximations involved in this model 
but we believe that the result are more acurate than previousely obtained. 

The phase boundaries for various anisotropics were obtained. It was shown that for finite 
anisotropy, when Josephson coupling plays a role, it tends to smoothen out energy jumps 
compared to the case when only the electromagnetic coupling is kept. Thus energy jumps 
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FIG. 12: Comparison between two models: (a) Phase boundaries for two different models, this 



paper (squares), Ref. 



25l (circles) for anisotropy 7 = 125 (solid squares and circles) and 7 = 250 



(blank squares and circles), (b) ^(Qi) vs. temperature for 7 = 150. this paper (solid squares), Ref. 



25|(open circles) (c) and (d) show E and N^/N with the same symbols and 7 value as in part (b) 



tend to increase with increasing anisotropy. They are also found to decrease with increasing 
magnetic field along the phase boundary. The phase diagram clearly shows a shift in the 
transition temperature as a function of the anisotropy. 

For finite anisotropy up to a value of 500 we showed that the data can be collapsed to 
a straight line in the In-ln plot of i?7^ versus kT/eod meaning that the transition occurs 
when a single variable combination of temperature, magnetic field and anisotropy becomes 
critical. From here we could deduce the value of the Lindemann parameter to be about 0.3. 
For infinite anisotropy we obtained scaling when plotting vs. kT/eod. 

While keeping just the electromagnetic coupling we observe a typical 2D melting behav- 
ior towards high magnetic fields where the transition temperature T2D ~ 16 K becomes 
independent of the magnetic field. Our results for 7 = 00 were compared with a recent 
simulation study using the substrate model^^*^. The agreement between the two results is 
excellent. 
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Finally we have c omp ared the two cases where the electromagnetic interaction is treated 
approximately (Ref. 123 ) and exactly (this paper) for a value of 7 = 125. It is shown that in 
this case transition temperatures are very close and keeping just the in-plane interaction in 
the form of a modified Bessel function of the second kind is a good approximation, justifying 
the claim that we have made in a previous paper^^^. In that paper our main goal was to 
extract the effect of columnar pins on the melting transition. For larger anisotropics the 
deviations of the two treatments become more pronounced. 

Comparing with experimental results the simulations seem to agree better with the posi- 
tion of the so-called irreversibility line than with the position of the true melting line. This 
may be due to the effect of point defects in the experimental samples. Even pristine sample 
contains point defects that tend to shift the phase boundary towards lower temperatures 
and flatten it out somewhat at the low temperature side. Point defects can also have an 
effect on the energy and entropy jumps across the phase boundary. If one would like to 
extract the value of the anisotropy of the experimental sample from the coparison with the 
present simulations, we would estimate it to be in the 250-450 range. 
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APPENDIX A: STRAIGHT FLUX-LINE APPROXIMATION 



In this section we show that the interaction of a pancake with a stack of pancakes located 
at a distance R away from the given pancake and hned up along the z direction reduces to 
the modified Bessel function of the second kind, Kq{R/\). 

Consider the interaction of a pancake situated at z=0 with a stack of pancakes. We begin 
with Eq. ()2.20j) and ()2.21|) . The total interaction of the pancake with the FL will be given 
by 



Utotai [^]=2deoln^ + YJ^]+ ( ^) , (Al) 



XJ ^ R \X J \X 

where Yi{R/X) and Y2{R/X) are given by 



Yi 



(Peq 



exp(— |m(i|/A) j In 



R' 



(A2) 



and, 



m=—oo ' ^ 



R 



We have a geometric series in Eq. ()A2|) which can easily be summed over to give 

(R\__feo ( l + eM-d/X) \ C 
'[xJ- X \l-eM-d/X))^''R- 

As A/c? ~ 120 so we can approximate Eq. ()A4|) with 



^'(f) =-2^^oln^,. (A5) 

Replacing the summation over m to an integration over z and changing the order of inte- 
gration we can express Y2 {R/ X) as 
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With a simple change of variable over z, z = R' t and then again exchanging the order of 
integration we get 



(A7) 



Finally integrating over R' one obtains 



= 2rf£oi^o(y). (A9) 



Combining Eq. (IAT|) . (IA5jl and (IA9jl we get 



a; \\ 

APPENDIX B: ENERGY SUM OVER THE IMAGES 

We again consider a rhombically shaped region with side L and angle 9, unit vectors 
are ei and 62 with ei • e2 = cos 9, as was done in the appendix of a previous paper.— The 
Green's function Gq which describe the 2D coulomb interaction between one vortex and 
another including all its images, as is implied by the periodic boundary conditions is given 
by the solution to London's equation (see e.g. Ref. 

(1 - A2V^)Go(R, A) = 27rA25(R), (Bl) 
with the parameter A setting the scale for the ran ge o f the interaction. The solution of the 



above equation was derived in the appendix of Ref, 



25. The result was 



r m \\ - '^^^^ cos(t2n - 2Ti(3n) sinh(7„ti) + cos{t2n) sinh(7„(27r - tj)) 

C.olK,Aj- 2 2^ 7„(cosh(27r7„) - cos(27r/3„)) ' ^""'^ 



n=— 00 



where 



27ri?i 2n 



ti = — t2 = —{Ricos9 + R2), Pn = ncos9, 7„ = sin^ + LV(27rA)2. (B3) 

1j 1j 
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To get a formula for Logarithmic interaction all we have to do is to take the limit A tends to 
infinity. In this limit there is only one term which diverges in the series given above namely 
that corresponding to n = 0. Separating out the term corresponding to n = we get 



+00 



Go(R, A) = termo + sin 9 
where, 



n=l 



cos(t2^ — 27r/9„) sinh(7„ti) + cos{t2n) sinh(7„(27r — ti)) 
7„(cosh(27r7„) - cos(27r/5„)) 



term. 



sin 9 sinh(7oi(:i) + sinh(7o(27r — ti)) 







2 7o(cosh(27r7o) - 1) 

In the equation above 'Jq = L/ (2,71 X) tends to zero. Taking the limit 70 



(B4) 



(B5) 



gives 



ternit 



sin 9 







(B6) 



Now the first term in the equation above is a constant (albeit infinite) independent of Ri 
and R2 and after dropping it we are left with the following expression: 



term, 



sin 9 







-TT 



1-6 



R^ 
L 



6 



L 



(B7) 



Eq. ()B4|) with with 7„ = sin 9 together with Eq. ()B7|) constitute the required summation 
over a logarithmic potential under 2D periodic boundary conditions. 
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